In [23]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import nengo
from nengo.utils.matplotlib import rasterplot
import seaborn as sns
sns.set_style("white")
sns.set_style("ticks")
plt.rcParams['figure.figsize'] = (9.0, 5.0)
plt.rcParams['lines.linewidth'] = 2
plt.rcParams['xtick.labelsize'] = 'large'
plt.rcParams['ytick.labelsize'] = 'large'
plt.rcParams['axes.titlesize'] = 'x-large'
plt.rcParams['axes.labelsize'] = 'x-large'
def prettify():
    sns.despine()
    plt.tight_layout()
def video(fname, mimetype):
    from IPython.display import HTML
    video_encoded = open(fname, "rb").read().encode("base64")
    video_tag = '<video controls alt="test" src="data:video/{0};base64,{1}">'.format(mimetype, video_encoded)
    return HTML(data=video_tag)

How to build a brain
with Python

Trevor Bekolay

www.bekolay.org/pycon2015

$\quad\quad\quad\quad$ $\quad$ $\quad$

Neural simulators

  1. Emulate neurobiology.
  2. Operate in continuous time.
  3. Communicate with spikes.

In [85]:
!neurondemo
NEURON -- VERSION 7.3 ansi (1078:2b0c984183df) 2014-04-04
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2014
See http://www.neuron.yale.edu/neuron/credits

loading membrane mechanisms from /usr/local/Cellar/neuron/7.3/share/nrn/demo/release/x86_64/.libs/libnrnmech.so
Additional mechanisms from files
 cabpump.mod cachan1.mod camchan.mod capump.mod invlfire.mod khhchan.mod mcna.mod nacaex.mod nachan.mod release.mod
oc>
In [2]:
import neuron

soma = neuron.h.Section()
soma.L = 40  # Length in um
soma.diam = 20  # Diameter in um
soma.insert('hh')

dendrite = neuron.h.Section()
dendrite.L = 150
dendrite.diam = 3
dendrite.insert('pas')

dendrite.connect(soma, 0, 1);
In [3]:
stim = neuron.h.IClamp(soma(0.5))
stim.delay = 0
stim.amp = 1  # Amplitude in nA
stim.dur = 500

rec_t = neuron.h.Vector()
rec_t.record(neuron.h._ref_t)

rec_v = neuron.h.Vector()
rec_v.record(soma(0.5)._ref_v);
In [4]:
neuron.h.finitialize(-60)
neuron.init()
neuron.run(500)
In [5]:
times = np.array(rec_t)
voltages = np.array(rec_v)
In [6]:
plt.plot(times, voltages)
plt.xlabel("Time (ms)")
plt.ylabel("Voltage (mV)")
plt.xlim(right=500)
prettify()

TITLE Anomalous rectifier
UNITS {
        (mA) = (milliamp)
        (mV) = (millivolt)
} 
NEURON {
        SUFFIX Khcvode
        USEION k WRITE ik
        RANGE  gkbar, gk, ik
} 
INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
PARAMETER {
        v (mV)
        celsius = 37 (degC)
        mon = 1
      man = 1
      nan = 1
        gkbar = .0003 (mho/cm2)
        ekcvode    = -30 (mV)
}
DERIVATIVE state {  :Computes rate and other constants at current v.
                    :Call once from HOC to initialize inf at resting v.
        q10 = 3^((celsius - 37)/10)            
                :"m" potassium activation system
        minf = 1/(1+exp((v+78)/7))
      mtau= 38/q10
        ntau = 319/q10
        m += mon * ((man * 0.8 * (minf-m)/mtau) + (nan * 0.2 * (minf-m)/ntau))
}

In [7]:
import brian2 as br

# Parameters
C = 281 * br.pF
gL = 30 * br.nS
taum = C / gL
reset_v = -70.6 * br.mV
VT = -50.4 * br.mV
DeltaT = 2 * br.mV
Vcut = VT + 5 * DeltaT
tauw = 144*br.ms
a = 4*br.nS
b = 0.0805*br.nA
Vr = -70.6*br.mV

eqs = """
dvm/dt = (gL*(reset_v - vm) + gL*DeltaT*exp((vm - VT)/DeltaT) + I - w)/C : volt
dw/dt = (a*(vm - reset_v) - w)/tauw : amp
I : amp
"""
In [8]:
neuron = br.NeuronGroup(
    1, model=eqs, threshold='vm>Vcut', reset="vm=Vr; w+=b")
neuron.vm = reset_v
neuron.I = 1.5*br.nA
In [9]:
voltage = br.StateMonitor(neuron, 'vm', record=0)
spikes = br.SpikeMonitor(neuron)
br.run(500 * br.ms)
In [10]:
vm = voltage[0].vm[:]
for t in spikes.t:
    i = int(t / br.defaultclock.dt)
    vm[i] = 20*br.mV
In [11]:
plt.plot(voltage.t / br.ms, vm / br.mV)
plt.xlabel('Time (ms)')
plt.ylabel('Voltage (mV)')
prettify()

Scaling from neurons to brains

Option 1: random

Option 2: statistics

Option 3: function

$\quad$

In [12]:
with nengo.Network(seed=10) as net:
    ens = nengo.Ensemble(15, dimensions=1)
sim = nengo.Simulator(net)

activities, rates = nengo.utils.ensemble.tuning_curves(ens, sim)
plt.plot(activities[:, 0], rates[:, 0])
plt.ylabel("Activity of neuron (Hz)")
plt.xlabel("Pressure (arbitrary units)")
plt.ylim(top=400)
sns.despine()
plt.savefig('img/neuralcode-1.svg')

Neural code $\approx$ binary code

\begin{align*} x &= 13 \\ a &= \mathtt{\LARGE{1 1 0 1}}\\ d &= 2^3 \, 2^2 \, 2^1 \, 2^0 \\ \hat{x} &= 1 \times 2^3 + 1 \times 2^2 + 0 \times 2^1 + 1 \times 2^0 \\ &= 8 + 4 + 1 = 13 \end{align*}
\begin{align*} x &= 0.5 \\ a &= 23 \text{Hz}, 60 \text{Hz}, 3 \text{Hz}, \text{...} \\ d &= d_0, d_1, d_2, ... \\ \hat{x} &= 23 \times d_0 + 60 \times d_1 + 3 \times d_2 + \text{...} \\ &\approx 0.5 \end{align*}

Nengo is a neural compiler


int x = 3;
int y = x * x;
printf("%d squared is %d", x, y);
$$\Downarrow$$

...
    subq    $16, %rsp
    leaq    L_.str(%rip), %rdi
    movl    $3, -4(%rbp)
    movl    -4(%rbp), %eax
    imull    -4(%rbp), %eax
    movl    %eax, -8(%rbp)
    movl    -4(%rbp), %esi
    movl    -8(%rbp), %edx
    movb    $0, %al
    callq    _printf
    movl    $0, %edx
    movl    %eax, -12(%rbp)
...

val = Node(output=0.5)
x = Ensemble(n_neurons=40, dimensions=1)
squared = Ensemble(n_neurons=40, dimensions=1)
Connection(val, x)
Connection(x, squared, function=lambda x: x * x)
$$\quad$$ $$\Downarrow$$ $$\quad$$

Ensemble 'x'
------------
neuron_type=LIF()
gain=[   11.25053668   202.13530672     7.22401682    33.2550193     35.97808839
     8.09500913    18.48783009     4.52011101    31.35391507    22.79914949
    99.57712634     9.9972526     10.4300937     26.44722378    16.64511324
     8.56703477   204.64658142    20.72961175    12.59879579     4.14133187
    67.31407332    23.71409087     6.91468826    29.96178344    66.86728154
    11.62997747    65.68513384    26.42944824   515.40393534     9.28512135
     5.40755927    54.51278395    30.00517318  2038.11571883    12.07297056
     6.00363065    20.58459027    10.03664157    11.29825745    35.45346505]
bias=[  9.44466013e-01  -1.90908878e+02   7.95278618e+00  -1.54254303e+00
  -6.26779005e+00   8.23853682e+00  -2.61010180e+00   4.50550452e+00
   8.47744978e+00   6.02264271e+00  -7.36539327e+01  -1.48795118e+00
   2.60012838e+00  -1.00077690e+00  -4.52833925e+00  -7.33437033e-02
  -1.73268693e+02  -9.16093383e+00   1.29271729e+01   4.29083492e+00
  -5.35103933e+01  -5.36784048e+00   4.45757505e+00  -8.74177689e+00
  -4.05562843e+01  -2.98500272e+00  -3.48752843e+01  -1.68317917e+01
  -4.96956302e+02   1.02301275e+01   3.04630174e+00  -2.09019904e+01
  -8.92324954e+00  -2.01891216e+03   4.63307652e+00   6.01789071e+00
   8.42017799e+00   5.07842705e+00   8.50458248e+00  -1.77831709e+01]


Ensemble 'squared'
------------------
neuron_type=LIF()
gain=[ 166.98348239   13.55000393   11.10501859   53.70178142   24.2046276
   10.55205272   26.97733589   10.69131454  156.05276375   20.9357595
    5.43800614   68.57201592   72.65123458   33.9680321    24.57512569
   21.1976429     8.67434404   14.53680377    8.09275949   22.9356037
   16.37519448   22.30781058   45.03301059   33.3462742    51.18795506
    8.69509189   18.11261709   14.81518274   17.63978738   35.22776264
    7.22343981   13.65154979    9.23444351   23.98063544   49.24764976
   59.66650699   41.04237034   24.87989631  238.85308633    7.96430909]
bias=[-130.53655762    7.9815554     7.90733762  -46.13916876   -7.32546692
    8.29221574   -4.50890614   10.42385029 -130.8065894     2.41053476
    2.37360306  -54.36523527  -46.73389939  -22.48866651   -5.65097727
   12.08397309    9.58977002   -4.21937901   -0.69890639   -2.60386428
   14.00605036   -8.02315356  -17.02830542    2.50456163  -32.80185813
   -1.03692479   -0.87231805   10.27828142    8.38370787   -8.55839931
    4.0243432    -4.39155615    9.07749725  -10.17955921  -40.19829161
  -20.341822    -10.16668571    6.74012336 -226.96010992    6.54200717]


Connection 'x->squared'
-----------------------
decoders=[[  1.60794810e-04   1.28863058e-04  -1.58415261e-05   2.33935353e-04
    2.27420025e-04  -1.99271624e-05   2.04197749e-04   2.64319388e-05
    1.48746335e-04   6.63795669e-05   2.90423062e-04   1.49912439e-04
    1.19707234e-04   1.92758931e-04   1.70101357e-04   1.47749115e-04
    4.91023627e-04   4.16415977e-04  -5.59905588e-05  -8.46265664e-05
    2.47756109e-04   2.65295259e-04   5.60745015e-05   2.08307183e-04
    2.63151063e-04   1.43506596e-04   2.64467801e-04   1.86141162e-04
    1.37304126e-04  -1.54465394e-04   7.12132301e-05   2.35052901e-04
    3.17766505e-04   4.28723382e-05   6.92211091e-05  -1.28598294e-04
    3.09857283e-05   1.16231202e-04  -4.20618385e-05   5.46221120e-04]]
solver_info={'rmses': array([ 0.022518])}
In [13]:
import nengo

with nengo.Network(seed=14) as net:
    val = nengo.Node(output=0.5)
    x = nengo.Ensemble(40, dimensions=1)
    squared = nengo.Ensemble(40, dimensions=1)
    nengo.Connection(val, x)
    nengo.Connection(x, squared, function=lambda x: x * x)
    
    v_pr = nengo.Probe(x.neurons[0], 'voltage')
    sp_pr = nengo.Probe(x.neurons, 'spikes')
    x_probe = nengo.Probe(x, synapse=0.03)
    sq_probe = nengo.Probe(squared, synapse=0.03)
In [14]:
sim = nengo.Simulator(net)
sim.run(1.0)
Simulation finished in 0:00:01.                                                 
In [ ]:
print("Ensemble 'x'")
print("------------")
print("neuron_type=%r" % x.neuron_type)
print("gain=%s" % sim.data[x].gain)
print("bias=%s" % sim.data[x].bias)
print("\n")
print("Ensemble 'squared'")
print("------------")
print("neuron_type=%r" % squared.neuron_type)
print("gain=%s" % sim.data[squared].gain)
print("bias=%s" % sim.data[squared].bias)
print("\n")
print("Connection 'x->squared'")
print("-----------------------")
print("decoders=%s" % sim.data[net.connections[-1]].decoders)
print("solver_info=%s" % sim.data[net.connections[-1]].solver_info)
In [15]:
plt.plot(sim.trange(), sim.data[v_pr])
plt.ylabel("Voltage (arbitrary units)")
plt.xlabel("Time (s)")
prettify()
In [16]:
rasterplot(sim.trange(), sim.data[sp_pr])
plt.ylabel("Neuron")
plt.xlabel("Time (s)")
prettify()
In [17]:
plt.plot(sim.trange(), sim.data[x_probe], label="$x$")
plt.plot(sim.trange(), sim.data[sq_probe], label="$x^2$")
plt.ylabel("Decoded value")
plt.xlabel("Time (s)")
plt.legend(loc="best", fontsize="x-large")
prettify()
In [18]:
val.output = lambda t: np.sin(t * 2 * np.pi)
sim = nengo.Simulator(net)
sim.run(1.0)
Simulation finished in 0:00:01.                                                 
In [19]:
plt.plot(sim.trange(), sim.data[v_pr])
plt.ylabel("Voltage (arbitrary units)")
plt.xlabel("Time (s)")
prettify()
In [20]:
rasterplot(sim.trange(), sim.data[sp_pr])
plt.ylabel("Neuron")
plt.xlabel("Time (s)")
prettify()
In [21]:
rasterplot(sim.trange(), sim.data[sp_pr])
plt.ylabel("Neuron")
plt.xlabel("Time (s)")
prettify()
plt.savefig('img/raster.svg')
In [22]:
plt.plot(sim.trange(), sim.data[x_probe], label="$x$")
plt.plot(sim.trange(), sim.data[sq_probe], label="$x^2$")
plt.ylabel("Decoded value")
plt.xlabel("Time (s)")
plt.legend(loc="best", fontsize="x-large")
prettify()
In [26]:
video("img/spaun.mp4", "mp4")
Out[26]:
In [ ]:
bg = nengo.networks.BasalGanglia(dimensions=4)
In [ ]:
import nengo_ocl
ocl_sim = nengo_ocl.Simulator(network)

import nengo_spinnaker
spinnaker_sim = nengo_spinnaker.Simulator(network)
In [168]:
%%writefile oscillator.py
import nengo
import nengo_viz
import numpy as np

tau, r = 0.01, 4

def feedback(x):    
    return [-tau*r*x[1]+1.01*x[0], tau*r*x[0]+1.01*x[1]]
def osc_shape(x):
    theta = np.arctan2(x[1], x[0])
    r = (2 - 2 * np.sin(theta)
         + np.sin(theta)*np.sqrt(np.abs(np.cos(theta)))
         / (np.sin(theta)+1.4))
    return -r*np.cos(theta), r*np.sin(theta)

with nengo.Network('Oscillator', seed=1) as model:
    stim = nengo.Node(0)
    oscillator = nengo.Ensemble(1000, dimensions=2)
    shape = nengo.Ensemble(100, dimensions=2, radius=4)
    nengo.Connection(stim, oscillator[0])
    nengo.Connection(stim, oscillator[1])
    nengo.Connection(oscillator, oscillator,
                     function=feedback, synapse=tau)
    nengo.Connection(oscillator, shape,
                     function=osc_shape, synapse=tau)

if __name__ == '__main__':
    nengo_viz.Viz(__file__).start()
Overwriting oscillator.py

Benchmarks

In [57]:
x = np.arange(8)
bit2 = np.array((x % 8) >= 4, dtype=float)
bit1 = np.array((x % 4) >= 2, dtype=float)
bit0 = np.array((x % 2) >= 1, dtype=float)
In [69]:
plt.plot(x, bit0 - 0.02, drawstyle="steps-mid", label="$2^0$")
plt.ylabel('Activity of bit (True or False)')
plt.xlabel('Encoded value ($x$)')
plt.ylim(top=1.2)
plt.legend(loc='lower center', fontsize='large', ncol=3)
sns.despine()
plt.savefig('img/binarycode-1.svg')
In [70]:
plt.plot(0,0)
plt.plot(x, bit1 - 0.02, drawstyle="steps-mid", label="$2^1$")
plt.ylabel('Activity of bit (True or False)')
plt.xlabel('Encoded value ($x$)')
plt.ylim(top=1.2)
plt.legend(loc='lower center', fontsize='large', ncol=3)
sns.despine()
plt.savefig('img/binarycode-2.svg')
In [71]:
plt.plot(0,0)
plt.plot(0,0)
plt.plot(x, bit2 - 0.02, drawstyle="steps-mid", label="$2^2$")
plt.ylabel('Activity of bit (True or False)')
plt.xlabel('Encoded value ($x$)')
plt.ylim(top=1.2)
plt.legend(loc='lower center', fontsize='large', ncol=3)
sns.despine()
plt.savefig('img/binarycode-3.svg')
In [72]:
plt.plot(x, bit0 + 0.02, drawstyle="steps-mid", label="$2^0$")
plt.plot(x, bit1, drawstyle="steps-mid", label="$2^1$")
plt.plot(x, bit2 - 0.02, drawstyle="steps-mid", label="$2^2$")
plt.ylabel('Activity of bit (True or False)')
plt.xlabel('Encoded value ($x$)')
plt.legend(loc='lower center', fontsize='large', ncol=3)
sns.despine()
plt.savefig('img/binarycode.svg')
In [93]:
plt.plot(activities, rates)
plt.ylabel("Activity of neuron (Hz)")
plt.xlabel("Pressure (arbitrary units)")
sns.despine()
plt.savefig('img/neuralcode.svg')

Neural code $\approx$ binary code